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1. INTRODUCTION 


Over the last two decades, significant attention has been devoted to the 
development of lightweight, durable thermal protection systems (TPS) for 
future space transportation systems. Research programs are currently under- 
way at the Langley Research Center to investigate various metallic TPS 
concepts [1]. One of the proposed candidates is the titanium multiwall tile 
(see 12] and references therein for a discussion). Early design procedures of 
the TPS concept involved both analytical and experimental studies. In 
particular, a degree of confidence has been established in the TPS concept due 
to the design studies by Jackson and Dixon [3] and Blair et al. [4], 

A titanium multiwall tile consists of alternating layers of superplas- 
tically formed dimpled sheets and flat septum sheets of titanium foil. As de- 
scribed in reference [3], this multiwall concept impedes all three modes of 

heat transfei: conduction, radiation and convection. The superplastically 

formed dimpled sheets and the long thin conduction path tend to minimize heat 
conduction. The flat septum sheets of titanium foil impede radiation. The 
small individual volumes created by the dimpled layers virtually eliminate air 
convection. The optimal design of such thermal protection systems requires 
effective techiques in coupled thermal and stress analyses. Finite element 
methods offer the greatest potential in modeling such complicated problems. 
However, the resulting semi-discrete equations may Involve many thousand 
degrees of freedom. Since the problem to be solved is transient and non- 
linear, the selection of an appropriate time integration method is an essen- 
tial step in the solution of such a complicated problem. AdeLman and Hafka 
[5] recently conducted a survey study on the performance of explicit and 
implicit algorithms for transient thermal analysis of structures. Calcula- 
tions were carried out using the SPAR finite element computer program [6] and 
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a special purpose finite element program incorporating the GEARB and GEARIB 
algorithms. Based upon their studies, they concluded that, generally, implic- 
it algorithms are preferable to explicit algorithms for "stiff" problems, 
though non-convergence and/or wide-banding of the resulting matrix equations 
may decrease the advantage of the Implicit methods. 

These difficulties are similar to those found in fluid-structure prob- 
lems. Over the past few years, several remedies have been proposed for these 
difficulties. Belytschko and Mullen [7] have proposed an explicit-implicit 
method where the mesh is partitioned into domains by nodes and the partitions 
are simultaneously integrated by explicit and implicit methods. Hughes and 
Liu [8] have proposed an alternate implicit-explicit finite element method 
where the mesh is partitioned into domains by elements and this element parti- 
tion concept simplifies the computer-implementation and enhances its compati- 
bility with the general purpose finite element software. 

Although the implicit-explicit method has been proven to be very success- 
ful in some fluid-structure Interaction problems (see e.g. , [8-10]), the size 
and complexity of the program are increased because of the addition of the 
implicit method. To overcome these difficulties, Belytschko and Mullen [11] 
have proposed an e’’^-E partition, in which explicit time integration is used 
throughout. However, different time steps within different parts of the mesh 
can be employed simultaneously. Partitioned and adaptive algorithms for ex- 
plicit time integration have also been proposed by Belytschko [12]. 

Recently, Liu and Belytschko [13] put forward a general mixed time 
implicit-explicit partition procedure within a linear context. It incorpo- 
rates the mentioned algorithms as special cases and is shown to have better 
stability properties than that in E™-E partition [11] . Similar concepts can 
also be used in transient conduction forced-convection analysis (see Liu and 
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Lin [14] ). 

In the present report, these implicit-explicit concepts ( nodes and 
elements) are extended to transient thermal analysis of structures where 
different time integration methods with different time steps can be used in 
each element group. The aim of this approach is to achieve the attributes of 
the various time integration methods. 

For example, in transient structural analysis, explicit methods require 
the size of the time step to be proportional to the length of the shortest 
element; while in transient thermal analysis, explicit methods require the 
step size to be porportional to the the square of the length of the shortest 
element. So it is more advantageous to employ this mixed time implicit- 
explicit technique for transient thermal analysis of structures since the 
E™-E partition proposed in [11,12] is often inefficient for this kind of 
problem though it is very efficient in structural analysis. 

In section 2 the finite element formulation for transient heat conduction 
is reviewed. In section 3 the mixed time integration procedures viz two 
element groups "A" and "B" are described. A family of integration partitions 
can then be deduced by selecting the appropriate definitions for the 
quantities of "A" and ”B". Five useful partitions which are of practical 
importances are presented. The stability criterion and critical time step 
estimates are given in sections 4 and 5, respectively. In section 6 the mixed 
time methods described in section 3 are generalized to NUMEG element groups. 

A computational algorithm for this mixed time implicit-explicit integration is 
also presented. In section 7 an illustrative example problem is described to 
demonstrate the practicability and usefulness of the proposed approach. In 
addition, the selection of a time integration method and the selection of an 
element group time step for each group are illustrated. Three numerical 
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example solutions are presented in section 8 to evaluate the performance 
(i.e., accuracy and stability behavior, computer storage and solution time, 
etc.) of these mixed time finite element algorithms. This represents the 
first comprehensive study of the effectiveness of the proposed methods. 
Related discussion and conclusions are presented in section 9. 
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2. FINITE ELEMENT FORMULATION FOR TRANSIENT HEAT CONDUCTION 

Consider a body enclosed by surface F which consists of two parts; 

r and r . Tlie Cartesian coordinates of the body will be denoted by x. . 
g q 1 

The governing equation for transient heat conduction with constant 
coefficients is: 


0,. , = - e 

il a 


in 0 , i=l , . . . , NSD 


( 1 ) 


6 = g 


for X. on r and t > 0 
i g 


( 2 ) 


0, .n. + he = q 
1 i 


for Xj on r and t > 0 

i q 


(3) 


and 


9=9 for X. in Q and t = 0 . (4) 

0 1 

Here a comma designates a partial derivative with respect to x^; a super- 
script dot designates time (t) derivative; n^^ is the x^ component of the 
outward unit normal vector; NSD is the number of space dimensions; a is the 
thermal diffusivity (the ratio of thermal conductivity to specific heat times 
density); 9 is the temperature; h is the ratio of convective heat transfer 
coefficient to thermal conductivity; and g, the prescribed boundary 
temperature; q, the prescribed heat flux and 9^, the initial temperature are 
given. Repeated indices denote summations over the appropriate range. 

The variational or weak form of equations (l)-(3) with (4) as the initial 


condition is: 
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(0 v) + A(9 v) = (q v)^ 

q 

where v is the test function; and 


(9 v) = — 0 vdn 


a 


(5) 


( 6 ) 


A(9 v) = e,^v,^d.Q + /j. h0v dr , (7) 

q 

and 

(q v)j, = /j, qv dr . (8) 

q q 

The finite element equations are obtained by approximating the trial 
functions by shape functions (N^) so that 

NEQ 

V = Z N [x. jd (t) , (9) 

i=l ^ J ^ 

NUMNP 

g = Z N [x jg [x ,tj , (10) 

l=NEQfl J J 

and 

0 = V + g (11) 

Here NUMiJP is the total number of nodal points used in the finite element 
mesh; and NEQ is the number of trial functions used (for this particular case 
it is equal to the number of equations to be solved). 

The resulting semidiscrete equation for transient heat conduction is 



- 7 - 


then: 

s 

M9 + Ke = F , 

with initial condition 

9(0) = 0^ , 

~ ~o 

where the mass matrix, 

M = = (N^ N.J 


= / 


U a 


N^N.dfi 
i J 


( 12 ) 

(13) 

(14) 


represents the thermal energy stored in the body The conductivity matrix, 


K = Ik. I = A[N, N J = r N, ,N. , dJ2 + /„ hN.N, dP 
~ ij-‘ ^ i i,k j,k ■'r 1 j 


(15) 


represents the conductive transfer of energy within the body and the 
convective transfer of energy across the boundary, The heat load vector, 

F = [Fj = [q Njj, - N^Jgj^- A(N^ N^jg^ . (16) 

q 


represents the impressed temperature condition on surface F and the surface 

S 

convection on surface F^. M and K are symmetric and positive definite. 
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3. MIXED TIME PARTITION PROCEDURES 

In this section, mixed time integration methods are employed to solve 

equations (12) and (13). For the purpose of describing these mixed time 

integration techniques the mesh is subdivided into element groups A and B, 

each of which is to be integrated by a different method. Let n be the time 

step number; 0 , V and F be approximations to 9(t ), 9(t ) and F(t ) re- 

~n ~n ~n ~ n ~ n ~ n 

spectively. Let mAt and At be the time steps used for element group A and 
element group B respectively, where m is an integer and is greater or equal to 
1. A time step cycle (mAt) can then be defined by an increment of m substeps 
with a time step of At each, so that one time step cycle is defined by step n 
to step rrhn. The portions of the matrices obtained by assembling element 
group A and element group B are denoted by superscripts "A” and ”B”, respec- 
tively. Hence it follows that any global matrix is the sum of the two ma- 
A B A B 

trices, cf. M = M + M and K = K + K . Nodes associated with only element 

^ ^ ^ m 

group B are denoted by superscript "B”, whereas those which are in contact 

with at least one element of group A are denoted by superscript "A"; nodes 

which are connected to both group A and group B are designated by "C", so "C" 

is a subset of ”A". To simplify the presentation, we further denote those 

element matrices associated with at least one node C are denoted by super- 
C G B "R 

script "C", so M and K are subsets of M and K respectively. However, in 
actual computer Implementation this element group is not necessary. With 
these definitions, M^= M^+ M^ and K^+ K^. 

Oi* ^ ^ 

Similarly, all vectors are then partitioned accordingly into "A" and"B" 
parts, cf. 9 = (6^9®)'*', V = (V^V^)'^ and F = (F^®)"^. The superscript "T" 

denotes the transpose. The vector 9 is sometimes redefined by augmented 

*A *R *A A T *R Pi T 

matrices, 9=9 +9 where 0 =(90) and 9 =(09). Similar defini- 

A 

tions are used for V and F. Any nonzero terms in F obtained in a corapu- 
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tation of F are neglected; they are assumed to be zero. 

As an example, consider a one dimensional bar with two groups of material 
A and B depicted in figure 1. Group B material has higher conductivity than 
group A that a smaller time step (At) is required for group B whereas a larger 
time step (mZ!it) can be used for group A. This mesh consists of 8 nodes and 7 
elements. Then the set of nodes "A" will be 1,2, 3, 4, 5; the set of nodes "B" 
will be 6,7,8; and the set of nodes ”C" will be 5. 


ELEMENT 6 = 

0 

© 

< 

(D © 

1 

1 

j ® © © 

NODE L = 

1 2 

3 4 

1 

5 6 7 8 

1 


GROUP A (mAt) 

1 

1 GROUP B (At) 
1 
1 


Fig. 1. One dimensional mesh with two element groups 


Let M^’, and F® be the thermal energy stored in the eth element, the 
conductive transfer of energy in the eth element including the convective 
transfer of energy across the element boundary and the heat load contributions 
to the global arrays respectively, then 




M 



Z- K 
e=”l ~ 




e=5 









•» 
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and 


M^= Z, 

~ e=l ~ 


K^= Z, K’ 
~ e=l 


,e 


F = Z, F 
~ e=l ~ 


*B *B 

If we let Pj, be the ith component of the global assembled vector P 


then 


*B 


P [0,0,0,0,P5=0,P^,P7,Pgj . 


With these definitions, the mixed time partition is given as follows. 
Governing equation 
for j=0,m; 


R^*A B^*B 

MV , , + + K“e 1, = F . , 

— n+j ~ ~n+j ~n+j ~n+j 


(17) 


and 


for J=l, . . . ,m-l; 


+ K®0*®. + K^e*^. = F*®. 


(18) 


ft \ ^X j / ^x 

where 9 , . is a suitable extrapolator (and/or interpolator) of 0 (and/or 9 . ) 

~n ~n+m 

for x=A and B. In actual computation, equation (18) is implicitly included 
in equation (17); and for j=l,...,m-l no quantities of A are being solved. A 


family of integration partitions can then be deduced from equations (17) to 
(18) if M is assumed to be lumped. Some members which are of practical im- 
portances are shown in table 1. 
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Table 1 



Time Integration In 

Time Integration In 

Ext rapolator/ Interpolator 

Designation 

Element Group A 

Element Group B 

Node A 

Node B 

E-E 

explicit with At 

explicit with At 

9*^ 

*B 

0 

~n 

mE-E 

explicit with mAt 

explicit with At 

9*^ 

~n 

e*® 

~n 

mE-T 

implicit with mAt 

explicit with At 

~nfm 

e*® 

~n 

E-I 

implicit with At 

explicit with At 


9*® 

I- 1 

implicit with At 

implicit with At 

Sn+l 



For purposes of describing the computer implementation and stability 
analysis, the modified generalized trapezoidal rule will be used to carry out 
the time temporary discretization of equations (17) and (18) though other 
implicit integration methods can also be used. 

® Modified generalized trapezoidal rule 
for. j=l, . . .. ,m; 


+ (l-a)jAt / 
~n+j ~n ~n 


(19) 


for 1 < j < m define the set "C" only, 


9® . = 0®, . , + (l-a)At V®, . T , 
j “■ 1 j “'1 


( 20 ) 


9^ =9^ + am At , 

~n+ra ~n+m ~n+m 


( 21 ) 
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and 


B® 

~n+j 


0®. . + aAt V®, . 
~n+j ~n+j 


( 22 ) 


In the above equations, a is a free parameter which governs the stability and 
accuracy of the method. Some useful partitions which have been depicted in 
table 1 are illustrated below. 


Example 1 : E-E Partition 

In this case, m=l, 9^,, = for x»A and B. Equations (17) to (18) 

~n+l ~n+l 

reduce to: 


MV . . 
— n+1 


-^+1 


5^+l~ £n+l ’ 
0 + (l-a)At V 




(23) 


(24) 


and 

In+l'^ ^+1 * 

Equations (23) to (25) represent the predictor-corrector explicit 
algorithms with equation (24) as the predictor and equation (25) as the 
corrector. 

Example 2 : mE-E Partition 

In this case, m > 1, 0^, , = 0'^, , and 0^, . H 0^ , . . Equations (17) to (22) 

~n+j ~-n+j ~n+j ~n+j 


reduce to: 
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PREDICTOR PHASE: 


equation (19) 


(26) 


and 


equation (20) 


(27) 


GOVERNING EQUATIONS: 


equation (17) (28) 


and 


equation (18) 


CORRECTOR PH/iSE: 


equation (21) 


(29) 


(30) 


and 

equation (22) (31) 

Example 3 : mE-I partition 

* A A 

In this case, m > 1; in equation (17), 9^^ = 9^^ for element group 

^A ~A C 

A only, 9 . = 9 for the portion which is related to K for j*0 and m. 
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This is done automatically if element group A is defined to be the implicit 

element group and element group B is defined to be the explicit element group. 

gA ^ gA 1 < j < m; and 6®, . = 0®. . for 1 < j < m. Equations (17) to 

~n+j ~n+j ~n+j ~n+j 

(22) reduce to: 


PREDICTOR PHASE: 


equation (19) 


equation (20) 


GOVERNING EQUATIONS: 


for j=0,m; 


m , . + K^9*^. + K®0 . . = F . . , 


for j=l, • . • ,m"l; 


+ K®9*®. + K^9*^. 


F*® 

~n+j ’ 
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CORRECTOR PHASE: 


equation (21) 


equation (22) 


Example 4 : E-I partition 

This is a special case of example 3. 


^ s’In+i - s:n+i 


9^=0 + (l-a)At V 

~n+l ~n ~n 


9 1 1 “ 9 1 1 + V ^ 
~n+l ~n+l ~n+l 


(36) 

(37) 

Equations (17) to (22) reduce to: 

(38) 

(39) 

(40) 


Equations (38) to (40) represent the implicit-explicit algorithms developed by 
Hughes and Liu (see e.g., [8,10]) in which equation (39) is the predictor and 
equation (40) is the corrector. 


Example 5 : I-I partition 

In this case, hf=1, 6'^,, = 0^,, and 6®,, 5 9^ , i • Equation (17) to (22) 

~n+l ~n+l ~n+l ~n+l 

reduce to the usual implicit formulation and it is: 


MV' . + K0 = F 
~~n+l ~~rrH ~n-(-l 


9 , = 9 + (l-a)At V 

~n+l ~n ~n 


9 , 1 = 9 . , + aAt V , . 
~n+l -^+1 ~n+l 


(41) 


(42) 


( 43 ) 
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4. STABILITY CRITERION 

The stability characteristics of these mixed time partition algorithms 
can be deduced using an energy balance technique (see [8] for a discussion). 
The analysis is restricted to the case in which F=0 and all capacitance 
matrices are lumped. To simplify the subsequent writing, the following 
notations will be used. 


x_. I = X - X , 

- -'n+m ^ ~n+m ~n 


(44) 




(45) 


I X , . I = X , . , , - X , . , 

'-~n+3 •' ~n+j+l ~n+j 


(46) 


and 


for j=0, 1, . . . ,m-l 


<x > = + X J/2 . 

~n+j ^~nd-j+l ~n+j 


(47) 


Assume : 


<V*;^ >'^ <V*^ >/z*® z*® > 1/ra^ , where 

~n+m ~ ~n+m ~ ~ 


*B m-2 *B 
z = V + <V “.> ; 

~ ~n j=l '^+1 


(48) 


<V*®>^ K*" <V*®>/V*'^ V*^ > (1 +(m-l)a)^ ; 

~n ~ ~n ~n ~ ~n 


( 49 ) 


*B .T *B . /„*A ,*A. .2 . 

<V , K <V , .>/V K V > (1-a) for 
~n+j ~ ~n+j ~n ~ ^n 
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( 50 ) 


and let: 


S . V*''’' V** + H*" V*“ > 0 ; 

n ~n ~ ~n ~n ~ ~n 


(51) 


> > 0 ; 

n-hn -'n-hn ~ ~-n+m 


(52) 


6 * 


P® , = <V*®. K® <v*j > > 0 ; 

n+j ~n+j ~ ~n+j 


(53) 


the energy expression of these mixed time partition procedures can be shown to 
be : 

S ^ < S - 2mAt - 2At P®, . . (54) 

n+m n n-hn j=0 n-tj 


Here, K = K - K and the stability is governed by M and M provided 
a > 1/2. Let 


2^ = ^ - 1/2 j At , (55) 

and 

Wj = ^ 4- 'a - l/2jjAt ; (56) 


the definitions of 



and M ° 


for the five cases discussed in section 3 are: 



Example 1 : E-E partition 


M 


*^= 2i . and 2? 


( 57 ) 


Example 2 : mE-E partition 


*R R B 

M = 0 , and M = 0, 


(58) 


Example 3 : mE-I partition 


M*^= 4- , and M*®= 


Example 4 : E-I partition 


M 


*R 


= hJ + 2^ > 2? 


(59) 


(60) 


Example 5 ; I-I partition 

M*^= , and M*^= (61) 

*R 

These mixed time partition procedures are stable if a > 1/2 and M 
*B 

and M are both positive definite. A summary of the results is as follows 


Example 1 : E-E partition 
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Example 2: iiiE-E partition 


= qWC ^2 and SJ® . <2 

crit crit crit 


(63) 


Example 3 : mE-I partition 


< 2 , and fi® , ^ < 2 . 

crit crit 


(64) 


Example 4 : E--I partition 


or , < 2 , and if f2 . 

crit crit 


(65) 


Example 5 : I-I partition 


unconditionally stable. (66) 


In equations (62) to (65), is defined to be jAt where 

denotes a typical eigenvalue of the eigenproblem 


M^e + K^e = 0 


(67) 
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S. CRITICAL TIMK STKI’ RSTIMATKS FOR BILINFAR QUADRILATERAL ELEMENT 

For the one point quadrature element, there are four eigenvalues of which 
only two of them are nonzero. One of these zero eigenvalues can be removed if 
the stabilization stiffness is added to the one point quadrature stiffness. 
This third eigenvalue can be shown to be bounded between the other two nonzero 
eigenvalues. In particular, for the case of a rectangular element, these 
three eigenvalues are Identical to those three obtained by using a two by two 
quadrature element. Hence a computationally-useful method of estimating the 
critical time step for linear quadrllaterlcal element can easily be derived 
for the mixed time methods. 

Following the previous section. At . < 2/X is the critical time step 

crit max 

restriction and X is the maximum eigenvalue of the e^^ element subjected to 
max 

insulated boundary conditions. 



Fig. 2(a). Arbitrary four-node quadrilateral element; 
(b) rectangular element 


For an arbitrary four node quadrilateral as shown in Fig. 2a, Let 
X = (x. ,x^,x-,x, ), 


and 
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I = (ypy2»y3.y4^’ 

are the x and y coordinates of the four nodes. Define 


and 


IJ 


Xi - Xj, 


for I, J=l, ... ,4 


Let 


^IJ 

" ^1 ■ yj’ 


2 

2 

X = 

^24 + %2’ 


2 

2 

Y =■■ 

^31 ^3’ 


^ “ ■^^ 24 ^ 31 -^ 42 ^ 13 ^’ 

and A is the area of the quadrilateral. It can be shown that 


At it < 2 — 

a{(X+Y)±/[(X-Y)^+4Z^J} 

where a is the thermal diffusivity and y is between 0.95 and 1.0. A numerical 
study of the eigenvalues of the various shapes (for all practical purposes) 
two by two quadrature elements has been performed. It is found empirically 
that y can be picked to be 0.95 if the element is really skewed and y can be 
picked to be 1.0 if the element is rectangular. In the special case of a 
rectangular element (i.e., with y=1.0), 


At 


crit 


„2 

!L Z 

- -^1 
2a’ 2aJ 




where Z and Z are defined in Fig. 2b. 
X y ® 


( 69 ) 
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6. IMPLEMENTATION ASPECTS 

In this section, the mixed time integration methods described in section 
3 are generalized to NUMEG element groups. Different time integration methods 
(implicit/explicit) with different time steps can be used in each element 
group. Let and be the element group time step and element group 

time respectively for NEG = 1,..., NUMEG. There are NUMEL elements in each 
element group, and At denote the minimum time step amount all these element 
groups. In this formulation all element group time steps are required to be 
integer multiples of At and the time steps for adjacent groups are integer 
multiples of each other. If a physical situation occurs dictating the use of 
adjacent implicit groups they must be combined into one group with the time 
step equal to the smaller of the two groups. In addition, for each implicit 
group, that element group time step must be greater than those of the adjacent 
explicit groups. The main advantage of this m^^ implicit - m 2 explicit - m^ 
implicit -... etc., technique is to minimize the semi -bandwidth of complicated 
problems especially in the three-dimensional case. To illustrate the idea, 
consider the one dimensional mesh shown in Fig. 3. It consists of NUMEG 
element groups and NUMNP nodes. In this case NUMEG is equal to 4 and NUMNP is 
equal to 12. Node 1 is assumed to be an essential boundary 


I( 

6At) 

i 1 

1 1 

1 1 

J E(2At) J 

f. fl 

1(4 At) 

1 

1 

1 

1 E(At) 

NODE 1 2 

3 

V r •/ 

4 5 6 7 

1 1 

8 9 

10 II 12 
1 

GROUP 

1 

1 1 

1 2 1 

1 1 

3 

1 4 

1 


Fig. 3. line dimensional mesh with four element groups 


condition node so that the number of equations, NEQ, is equal to 11. The 
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essence of the present development can be deduced graphically by considering 
the solution procedures of the matrix equations. The "active column equation 
solver" is the key to the success of this technique (see [8,13] for a 
description of this equation solver). The profile of the effective stiffness 

■k 

matrix K of this one dimensional mesh is shown in Fig. 4. It can be observed 


NODE NO 
2 

EOT NO 
1 

1 

2 

1 





T 


1 

3 

2 


3 

4 



GROUP 1 

(6A-t) I 

4 

5 

3 

4 



T 
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6 

5 

SYM 



I] 

GROUP 2(2At) E 

_L 1 

7 

6 



T 

n 

[i 

9] 





8 

7 







10 

TT] 




9 

8 

GROUP 3 (4A-t) I 

n 

12 

13 



10 

9 



_L 





n 

14 



II 

10 



T 







15 




KiKOUf ^ ) t 






12 

II 

1 


jL 







I6 


Fig. 4. Profile of the effective stiffness matrix 


from Fig. 4 the following: 


Group 1: implicit with 4t^= 6At , five words of storage (1-5), 3 elements and 
3 equations. 

Group 2: explicit with At^= 2At , two words of storage (6-7), 3 elements and 2 
equations . 

Group 3: implicit with At^= 4At , seven words of storage (8-14), 3 elements 
and. 4 equations. 
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Group 4: explicit with At^= At , two words of storage (15-16), 2 elements and 
2 equations. 

The equation systems of each element group are uncoupled and hence each group 
can be integrated at its own group time step. For example, assume the 
effective stiffness matrix K is formed and factorized once. Since each group 
has its own time clock that in a time interval of 6At, group 1 will be 
integrated implicitly once, group 2 will be integrated explicitly three times, 
group 3 will be Integrated implicitly once and group 4 will be integrated 
explicitly six times. In order to handle the forward reduction and 
backsubstitution and update procedures automatically, we required two arrays 
At„„_ and each has a dimension of NUMNP. At„_„„ array contains the 

nodal time steps of each node. Nodes associated with only one element group 
NEG are assigned a time step of whereas those which are in common to 

other element groups are assigned to have the maximum time step amount the 
adjacent groups. array contains the nodal time of each node. From 

these two arrays ^NODE^ boundary condition codes, another 

time step array and an equation time array of the equation systems 

can then be generated. Both At^^^ and have dimensions of NEQ. It is 

required further a master time which is incremented by the smallest time 
step At. For this particular example the Atj^^j^g, and At^^^^ arrays are: 

AtNODg“[ , 6At , 6At , 6At , 2At , 2At , 4At ,4At ,4At ,4At ,At ,At] ; 

and 

Atjj^Q=(6At,6At,6At,2At,2At,4At,4At,4At,4At,At,AtJ . 

The T,,_„„ and T.,„„ arrays are incremented by time steps of At.,„_ and At,,_^ 

respectively. With these definitions, the generalized mixed time integration 

is to proceed over the time interval [0,T ]. The procedures are as follows: 

max 
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1. Initialization. Set Tj^, Tj^gg, "^NODE '^NEQ”®’ define the 

initial data 0 and V 
~o ~o 

* 

2. Form and factorize K where 

* numeg 
K = Z K 
NEO 1 


NEMEL 

K = Z K 
e=l 


K ^ + aAt if Implicit, 

~ ~ NEG~ 


and 


K if explicit (lumped capitance matrix is assumed). 


3. Time step loop. 

^ ’’‘'max 

s A.- ^DISCRETE 

5. Set F = oiAt F 

NEQ ^ 

6. Loop on element groups NEG=i, . . . ,NUMEG; 

^*^NEG ■ 

7. Loop on elements e=l , . . . ,NUMEL; 

7a. Define predictor values 6®; 


r® 
NODE 


.e 

"NODE 


If T.. + AC “ < Torr, 


1 SuODE * “^'''nODE 'NODE' 

Otherwise f - + a-a)(T„-T=^^^)V»„^^. 


t 

means is replaced by", 
t t 

we use Torr=i.0E-9. 
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7b. 

and 

where 

7c. 

7d. 

6a. 

6b. 

8 . 

8a. 

8b. 

9. 

9a. 


9b. 

10 . 


^6 

From element effective force f 
f*e ^ implicit, 




f - = O" “ “ H explicit, 


W = diagonal matrix with At„ along the diagonals. 

Sum up effective force from elemental contributions 

* * 

F + F + f 

End of element loop . 

Update TjjgQ. 

T = T + At 

NEC NEC NEC 

End of element group loop. 

Solve for 9. (In order to enhance the efficiency of these mixed time 
methods, the active column equation solver used in [8,13] is modified 
in the forward reduction and backsubstitution procedures as follows;) 
Loop on equation number N=1,.,.,NEQ; 

“ “bEQ - 

Forward reduction and backsubstitution for equation N. 

End of equation number loop. 

Update V and £; 

Loop on N=l, . . . ,NUMNP; 

“ + ‘*SoDE - ■'« > S” 

9^ ■*- solution from step 8. 

V - (9 -9 )/oAt^Qj^g. 

End of nodal number loop. 

Update TjjQpg; 
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10a. Loop on N=1 , . . . ,MUMNP; 

If TnODE ^^NODE “ ^o 10b. 

e ^ 

NODE NODE NODE' 

10b. End of nodal number loop. 

11. Update Tjjgp; 

lla. Loop on N=1,...,NE0; 

1^ . t"" + At"" 

%EQ ^NEQ NEQ 

llb. End of equation number loop. 

3a. End of time step loop. 

A corresponding flowchart of the above procedures is shown on the next few 


pages 
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7. ILLUSTRATION OF THE SELECTIONS OF A TIME INTEGRATION METHOD AND AN ELEMENT 
GROUP TIME STEP 

A two dimensional finite element model for the Thermal Protection System 
(TPS) of the space shuttle [1] is formulated herein to illustrate the 
selection procedures for a time integration method and an element group time 
step. The practicability and usefulness of the algorithms can also be 
demonstrated with this model as will be described in section 8. Since the 
emphasis in this report is on the evaluation of the methodologies described in 
previous sections, nonlinear effects such as internal and external radiation 
are ignored. Further only the assumed mean temperature material properties of 
the various components of the TPS are considered. The finite element mesh is 
depicted in Fig. 5. Due to symmetry, only half of the TPS is modeled. Its 
material properties are tabulated in table 2. The number of elements, the 
minimum characteristic length !L . , the estimated explicit critical time step 
^^min’ proposed integration method and the proposed element group time 

step for each group are included in table 3. As can be seen from table 3 due 
to the various thermal time scales (e.g. At ^ =0.041 sec for AL, and 
At^n“l^*^’ sec for RSI), a single integration method is definitely not 
effective. For example, if an explicit method is employed, a time step of 
0.041 sec has to be used; while if an implicit method is employed, there is no 
stabilit 5 r-imposed limitation on At; however, wide-banding bandwidth and/or 
demanding computer storage of the resulting matrix equations proportional to 
the square of bandwidth may decrease its advantage. The family of mixed time 
integration schemes developed is best suited for this type of problem. The 
attributes of the various time integration methods are fully achieved using 
the proposed approach as can readily be seen from table 3. It should also be 
observed that At . can be set as high as 16.71 sec in this mixed model, the 
smallest time step required by an explicit group. Subsequently At=16.71 sec 
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is employed for element groups 1 and 2 (042 Coating and RSI), At=66.84 sec for 
element group 3 (RSI), and At=66.84 sec for element groups 4 and 5 (RTV, FKLT, 
AL and AIR). 

Hie advantage of this proposed ntLxed time implicit-explicit finite 
element concepts can further be visualized if a nonlinear analysis of the 
above finite element model is assumed. The thermal responses of the various 
components of the TPS may be divided into regions of slowly and rapidly 
varying temperatures and they are: 

1. 042 Coating 

TVie thermal dlffusivity is almost independent of temperature and the 
estimate At is so small (0.1225 sec) that explicit calculation is not 
cost effective at all. Implicit calculation is best suited since only 
limited partial reformation and refactorization of this element group's 
effective stiffness will be needed. This is one of the major advantages 
of the mixed time implicit-explicit concept. 

2. RSI 

The thermal dlffusivity changes rapidly with temperature, therefore 
reformation and refactorization of this element group's effective 
stiffness are frequently required if an implicit method is employed. 
Non-convergence and/or wide-banding of the resulting element group 
effective stiffness may decrease the advantage of applying the implicit 
method. The estimated At (based on mean temperature value) is 16.71 sec; 
therefore an explicit method is best suited (no matrices operations). 

3. RT7 

Its thermal dlffusivity is independent of temperature and the esti- 
mated At is so small (0.106 sec) that an implicit method is recommended. 



X MATERIAL 

X 

N 

X 

's 

PROPERTIES \ 

042 

COATING 

RSI 

RTV 

FELT 

AL 

AIR 

\ w/m-K 
conductivity 

1.4338 

0.1354 

0.3113 

0.0363 

13.5373 

0.0271 

p kg/m^ 
density 

1665.92 

144.17 

1409.62 

96.11 

2851.28 

1.126 

C J/kg-K 
specific heat 

1317.96 

1255.20 

1213.36 

1213.26 

962.32 

1009.0 

X 2, 

a= — m /sec 
pc 

thermal 
dif fusivity 

6.52xl0~^ 

7.48x10"^ 

1.82x10'^ 

3.11x10”^ 

4.78x10“^ 

2.39x10“^ 

The mean temperature material properties are computed from the 
1200 K, 950 K, 500 K, 477 K, 333 K and 300 K. 

average of those at 


































Table 3. Characteristic Length and Time Scales 


MATERIAL 

TYPE 

NUMBER OF 
ELEMENTS 


ESTIMATED 

EXPLICIT 

At 

(sec) 

042 

12 

0.04 

0.1225 

coating 




RSI 

30 

0.5 

16.71 

RSI 

18 

1.04 

72.29 

RTV 

18 

0.02 

0.106 

FELT 

6 

0.40 

25.70 

AL 

16 

0.20 

0.041 

i 

AIR 

20 

2.0 

8.36 


INTEGRATION ELEMENT GROUP 
METHOD TIME STEP 


total number of elements = 120 


I = implicit integration 
E = explicit integration 


_v 


16.71 


66.84 


66.84 
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4. 

Its thermal diffusivity is fairly independent of temperature (at least in 
the operational temperature range of the space shuttle) and its estimated 
At is so stringent (0.041 sec) that an implicit method is again 
recommended. 

5. FELT 

Its thermal diffusivity is independent of temperature and its estimated 
At is large (25.7 sec), therefore Implicit method or explicit method can 
be used. For this example, the implicit method is proposed. However for 
3D calculation, in order to reduce the bandwidth and the unnecessary 
nonlinear calculations (due to the fact that the adjacent groups might be 
implicit groups too), explicit method will be recommended. 

6. AIR 

The variations of the thermal diffusivity with temperature are small in 
our range of interest. Also the estimated At is small, therefore 
implicit method is again well suited. 

Remarks : 

(1) In the above At calculations, i . = min{£ ,£ } where l and Z are 

mxn X y x y 

defined in section 5. 

2 

(2) The estimated At is defined to be f , /2a. It is a safe critical time 

rain 

step calculation. See [13-14] for a discussion. 

(3) By virtue of the fact that each element group's effective stiffness is 
uncoupled to the global assembled matrix equations system, any element 
group can be reformed and refactorlzed at any instant if required without 
affecting the global equations system. This partial factorization 
procedure can further be enhanced if it is to be combined with an 
iterative update procedure. 
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8. NUMERICAL RESULTS 

A two dimensional finite element pilot computer code incorporating the 
methodologies described in previous sections has been written to evaluate the 
performance of these mixed time finite element algorithms. Three numerical 
examples are presented to demonstrate the^ accuracy, stability and efficiency 
of these proposed methods. All computations are performed on a CDC Cyber 
170/730 computer in single precision (60 bits per floating point word) and 
lumped capacitance matrices are used throughout. 

Example 1 Verification of the Modified Equation Solver 

The purpose of this numerical example is to confirm the generality of the 
modified column equation solver. The finite element mesh as shown in Fig. 6 
consists of 100 uniform square elements (each with £^=il^=10 m). The mesh is 
deliberately divided into four element groups; 

1) Group 1 is integrated implicitly with a group time step of 2At. 

2) Group 2 is integrated explicitly with a group time step of 2At. 

3) Group 3 is integrated implicitly with a group time step of 4At. 

4) Group 4 is integrated explicitly with a group time step of At. 



Fig. 6. Problem statement and finite element 
mesh of numerical example 1 
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The above partition is called a 2I-2E-4I-1E partition. Observed that 
group 3 (41) is chosen to be at the center of the square. As described in 
[8,13-14], the variable column heights (which affect the solution time) depend 
on the nodal numbering. Furthermore, no forward reduction and backsubstitu- 
tion are performed until every 4At time interval. As can be pictured from 
Fig. 6, the column heights of the "4At” equations are coupled to the "At" and 
"2At" equations, therefore unlike the one dimensional case (see [14] for a 
discussion) a special forward reduction and backsubstitution procedure as 
described in section 6 is required for efficiency. 

r\ 

The thermal diffusivity a is chosen to be 6.25 m^/sec for all the 
elements. Node 1 is kept at a constant temperature of 0.0 °C while all the 
other nodes are subjected to a constant initial temperature of 0.1 °C. All 
four sides are insulated. A time step of At equal to 1 sec is employed for 
the transient analysis. The temperature-time histories of nodes 3, 46, 57 and 
67 are depicted in Fig. 7. In order to confirm the accuracy of these 
solutions, the mixed time finite element solutions are compared to those 
obtained using a finite difference 21x21 (400 zones) grid with a single time 
step of At equal to 1 sec. The finite difference solutions of these four 
nodes are also depicted in Fig. 7. As can be seen, the comparison is quite 
good despite the crudeness of the mesh (as compared to the finite diference 
method) and the larger time steps used in the finite element mesh (1 as 
compared to 1, 2, and 4 for groups 4, 1 and 2, and 3 respectively). The 
maximum relatively difference between the two solutions is no more than 2%. 
Example 2 D emonstration of Accuracy, Stability and Efficiency 

The purpose of this numerical example is to evaluate the performance of 
this mixed time implicit-explicit method as compared to an implicit method. 

The problem statement is depicted in Fig. 8. The finite element mesh consists 
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of 400 uniform size rectangular elements (each has m and ^^“10 ®) 486 

nodes (81 in x direction and 6 in y direction) where the x-axis is defined by 
joining node 3 to node 17. The material properties of this finite element 
model are composed of four different thermal dif fusivities of 0.125a, 0.25a, a 
and 40a respectively as shown in Fig. 8. a is chosen to be 12.5 m /sec. In 
order to simulate a realistic three dimensional bandwidth using this two 
dimensional finite element model, the nodes are numbered sequentially along 
the x-axls. The heat load which is also shown in Fig. 8 is applied at node 
1. The initial temperature for all the nodes is 0.1 °C. All four sides are 
insulated. 

As mentioned earlier, due to the different thermal dif fusivities , a 
single time integration method is inefficient for this problem. For example, 
if explicit method is employed, a time step of 0.025 sec is required so that 
16,000 steps (as compared to 400 steps) are required for this problem. 
Therefore purely explicit integration method is not being considered here. If 
implicit method is employed, a much larger time step can be employed though 
the wide bandwidth and demanding computer storage (e.g. 3288 vs. 33771 single 
precision words) may decrease the advantage of the implicit method. This can 
readily be seen in table 4. Hence mixed time implicit-explicit method is best 
suited for this example. A total of five computer runs are made to 
demonstrate the accuracy, stability and efficiency of the proposed mixed time 
methods and they are: 

(1) 1E-4E-8I-8E with At=l sec. 

(2) II with At=l sec. 

(3) 1I-1E-2I-2E with At=4 sec. 

(4) II with At=4 sec. 


(5) II with At=8 sec. 
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Table 4. Storage (in Single Precision Words 


Item 

1E-4E-8E-8I 
(At=l sec) 

11 

(At=l sec) 

1I-1E-2E-2I 
(At=4 sec) 

11 

(At=4 sec) 

11 

(At=8 sec) 

Number of 
equations (NEQ) 

486 

486 

486 

486 

486 

Number of 
Terms in 
Stiffness (NA) 

3288 

33771 

6089 

33771 

33771 

Average Half 
Bandwidth 

H f* 

MB = NA/NEQ 

6 

69 

12 

69 

69 


Table 5. Solution Times (in Seconds) 
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In run 1, groups I, II and IV are Integrated explicitly with IE, 4E and 8E 
respectively; and group III is integrated implicitly. In run 3, groups IT and 
IV are integrated explicitly with IE and 2E respectively; while groups I and 
III are integrated implicitly with II and 2I respectively. The time 
step At and the element group time steps are chosen according to the critical 
time step formula given in section 5 with p equal to 1.0. The temperature- 
time histories at six different locations for run number 1, i.e. 1E-4E-8I-8E, 
are presented in Fig. 9. These results are virtually identical to those 
obtained from run number 2, i.e. II. It should be pointed out that the curves 
shown in Fig. 9 are plotted for every time interval of 4 except for run number 
5 in which the time interval is 8. Since the plot of node 327 is integrated 
with a time step of 8, therefore the temperature time history shown is "step- 
wise" (from time “ 48 to 240 sec). The storage requirements and solution 
times for the above five runs are tabulated in tables 4 and 5 respectively. 
From these tables, the advantages of these mixed time integration methods are 
apprarent. The relative accuracy of these five runs are compared and they are 
shown in Fig. 10. These curves are plotted at every time interval of 8 sec. 

In summary, from accuracy consideration, runs 1 and 3 made with the mixed time 
implicit-explicit method, are comparable to runs 2 and 4 made with an implicit 
method but with two different time steps. Run 5, also an implicit method but 
utilizing a larger time step, was not as accurate as the others. Runs 1 and 3 
required 10% and 18% respectively of the storage required by run 4 and 65% and 
50% respectively of computational time of run 4. The potential savings in 
both computer time and computer storage will be even greater in three 
dimensional and/or nonlinear calculations. Finally, the critical time step 
estimate is also confirmed by this numerical example. 






NODE 3 


-ie-4E-8I-8E 

AN01I(0T«1) 

- 1ME-2I-2E 
AN04I(0T*4) 

- 1I(DT = 8) 




TIME (SEC) 


Fig. 10. Comparison of sc 


TIME (SEC ) 


ions of five integration methods 
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Example 3 Applications to TPS 

The purpose of this numerical example has been described in section 7. 

Two integration methods are employed to compute the time histories of the top 
(node 84) and bottom (node 64) surface temperatures of the 042 coatings andthe 
lower (nodes 11 and 137) and upper (nodes 15 and 57) aluminum skin 
temperatures for the two dimensional finite element model described in section 
7. Its dimensions (not to scale and in cms) are depicted in Fig. 5. The heat 
loads which are also shown in Fig. 5 are applied at the top and bottom of the 
TPS model. The initial temperature for all the nodes is -20 °C. Sides AD and 
BC are symmetry planes. 

From the data given in section 7, 407561 time steps (as compared to 1000 
steps) are required if purely explicit method is employed. Therefore, purely 
explicit time integration method is not being considered. The 1I-1E-4E-4I-4I 
time integration with a time step of 16.71 sec for the 042 coating and part of 
RSI, and a time step of 66.84 sec for the other part of RSI, RTV, FELT, AL and 
AIR is proposed. The computed results are then compared to those using purely 
implicit method with a single time step of 16.71 sec. The temperature-time 
histories (plotted at every time interval of 66.84 sec) at six different 
locations (042 coatings and aluminum skin) are presented in Fig. 11. These 
results are virtually identical to those obtained from the purely implicit 
case, i.e. 11. The closed up temperature-time histories at four different 
locations are given in Fig. 12. These curves are plotted at every time 
interval of 16.71 sec for a time period of about 4000 sec. The solution times 
and computer storage requirements for these two runs are tabulated in tables 6 
and 7 respectively. Eventhough only a total of 147 nodes and 120 elements (a 
very small mesh as compared to a typical 3D finite element model of 3000 nodes 
and 4000 elements) is used, a factor of about 1,7 in solution time and a 
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factor of about 1.5 in computer storage are gained. Again, the advantages as 
well as the stability of these mixed time methods are fully illustrated. 














EMPERATURE 






Table 6. Solution Times (in Seconds) 


— IJetl^d 
Item 

1I-1E-4E-4I-4I 
(At=16.7l sec) 

11 

(At=16.71 sec) 

Formulation 
of Stiffness and 
Load 

105.12 

162.54 

Factorization 

0.36 

0.57 

Forward Reduction 
and 

Backsubstitutlon 

45.05 

91.00 

Total Solution 
Time 

160.42 

263.35 


Table 7. Storage (In Single Precision Words) 


-^ll^thod 

Item 

1I-1E-4E-4I-4I 
(At=16.71 sec) 

11 

(At=16.71 sec) 

Number of 
equations (NEQ) 

147 

147 

Number of 
Te nns in 
Stiffness (NA) 

2121 

2933 

Average Half 
Bandwidth 

1 f 

1 MB = nA/NEQ 

1 

14 

19 

1 
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9. CONCLUSIONS 

The computer Implementation aspects; stability criterion and the 
evaluation of the performance of mixed time implicit-explicit finite element 
as applied to multidimensional transient thermal analysis are presented. The 
report focuses on applications where the thermal model would normally have 
different time scales of thermal response. Three two-dimensional examples are 
presented to illustrate the approach: (1) comparison with finite difference 

method, (2) demonstration of accuracy, stability and efficiency and (3) 
application to Thermal Protection Systems. The examples show: (1) the 

superiority of mixed time methods to a single integration method (either 
implicit or explicit), (2) potential savings in computer time and computer 
storage, and (3) the accuracy and stability behavior of mixed time finite 


element. 



- 53 - 


REFERENCES 

1. H. N. Kelly; D. R. Ruramler; and R. L. Jackson: Research in Structures 

and Materials for Future Space Transportation Systems - An Overview. 
Presented at the AIAA Conference on Advanced Technology for Future Space 
Systems, May 8-11, 1979, Langley Research Center, Hampton, Virginia, AIAA 
Paper No. 79-0859. 

2. J. L. Shideler; H. N. Kelly; M. L. Blosser; and H. M. Adelman: Multiwall 

TPS - An Emerging Concept. AIAA/ASME/ASCE/AHS 22nd Structures, 

Structural Dynamics and Materials Conference, April 6-8, 1981, Atlanta, 
Georgia,. 

3. R. J. Jackson; and S. C. Dixon: A Design Assessment of Multiwall, 

Metallic Stand-Off and RSI Reusable Thermal Protection System Including 
Space Shuttle Applicatons. NASA TM 81780, April 1980. 

4. W. Blair; J. E. Meany; and H. A. Rosenthal: Design and Fabrication of 

Titanium Multiwall Thermal Protection System (TPS) Test Panels, NASA CR 
159241, February 1980. 

5. H. Adelman; and R. Haftka: On the Performance of Explicit and Implicit 

.Algorithms for Transient Thermal Analysis of Structures. NASA TM 81880, 
September 1980. 

6. SPAR Structural Analysis System Reference Manual, Vol. 1, NASA CR 158970- 
1, December 1978. 

7. T. Belytschko; and R. Mullen: Stability of Explicit-Implicit Mesh Parti- 

tion in Time Integration. International Journal Numerical Method in Eng . 
1^, 1575-1586, 1978. 

8. T. Hughes; and W. Liu: Implicit-Explicit Finite Elements in Transient 

toalysis. Journal Applied Mechanics 45, 371-378, 1978. 

9. W. Liu: Development of Finite Element Procedures for Fluid-Structure 

Interactions. EERL 80-06, California Institute of Technology, Pasadena, 
California, August 1980. 

10. T. Hughes; W. Liu; and A. Brooks: Review of Finite Element Analysis of 

Incompressible Viscous Flows by the Penalty Function Formulation. 

Journal Computational Physics 30 , 1-60, 1979. 

11. T. Belyt,schko; and R. Mullen: Explicit Integration of Structural 

Problems. In P. Bergan et. al. , (eds.) Finite Elements in Nonlinear 
Mechanic, 2, 697-720, 1977. 

T. Belytschko: Partitioned and Adaptive Algorithms for Explicit Time 

Integration. In the book of Nonlinear Finite Element Analysis in Struc- 
tural Mechanics, edited by W. Wunderlich et. al. , Springer-Verlag, 

Berlin, July 1980. 


12 . 



- 54 - 


13. W. Liu; and T. Belytschko: Mixed Time Implicit-Explicit Finite Element 

for Transient Analysis. Computers and Structures, Vol. 15, 1982, pp. 
445-450. 

14. W. Liu; and J. Lin: Stability of Mixed Time Integration Schemes for 

Transient Thermal Analysis. Numerical Heat Transfer Journal. Vol. 5, 
1982, pp. 211-222. 



1. Report No. 

NASA CR- 172209 


2. Government Acc^ion No. 


3. Recipient's Catalog No. 


4. Title and Subtitle 

MIXED TIME INTEGRATION METHODS FOR TRANSIENT THERMAL 
ANALYSIS OF STRUCTURES 


5. Report Date 

September 1983 


6. Performing Organization Code 


7. Author(s) 

Wing Kara Liu 


9. Performing Organization Name and Address 
Northwestern University 

Department of Mechanical and Nuclear Engineering 
Evanston, IL 60201 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546 


8. Performing Organization Report No. 


10. Work Unit No. 


11. Contract or Grant No. 
NAGl-210 


13. Type of Report and Period Covered 
Contractor Report 
10-1-81 - 9-30-82 


14. Sponsoring Agency Code 


15. Supplementary Notes 

Langley Technical Monitor: 
Final Report 


George C. Olsen 


16. Abstract 

The computational methods used to predict and optimize the thermal-structural 
behavior of aerospace vehicle structures are reviewed. In general, two classes 
of algorithms, implicit and explicit, are used in transient thermal analysis of 
structures. Each of these two methods has its own merits. Due to the different 
time scales of the mechanical and thermal responses, the selection of a time 
integration method can be a difficult yet critical factor in the efficient 
solution of such problems. 


Therefore mixed time Integration methods for transient thermal analysis of 
structures are being developed. The computer implementation aspects and numerical 
evaluation of these mixed time implicit-explicit algorithms in thermal analysis 
of structures are presented. A computatlonally-useful method of estimating the 
critical time step for linear quadrilateral element is also given. Numerical 
tests confirm the stability criterion and accuracy characteristics of the methods. 
The superiority of these mixed time methods to the fully implicit method or the 
fully explicit method is also demonstrated. 


17. Key Words (Suggested by Author(s)) 1 18. Distribution Statement 


Implicit Thermal response 

Explicit Critical time step Unclassified - Unlimited 

Implicit-explicit 

Mixed time Subject Category 64 

Transient thermal analysis 

19. Security Classif. (of this report) 20. Security Classif. (of this page) 21. No. of Pages 22. Price 

Unclassified Unclassified 55 A04 


N-305 


For sale by the National Technical Information Service, Springfield, Virginia 22161 



















